! 2002-9 (c) Mark A. Watson

MODULE fmm_main

   USE fmm_global_consts
   USE fmm_global_types

   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: fmm_initial,                      &
             fmm_final,                        &
             fmm_get_J_matrix,                 &
             fmm_get_multipole_potential

   ! The multipole potential
   REAL(REALK), POINTER, SAVE :: Vff(:,:)
   ! The (modified) multipole expansions
   TYPE(raw_mm_data), SAVE :: LHS_mms, RHS_mms

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE fmm_initial(mode,LMAX)

      USE fmm_scheme_module, ONLY: fmm_init_scheme
      IMPLICIT NONE
      CHARACTER(*) , INTENT(IN) :: mode
      INTEGER(INTK), INTENT(IN) :: LMAX

      INTEGER(INTK) :: job_type

      SELECT CASE (TRIM(mode))
         CASE ('MD4_FMM');  job_type = MD4_FMM
         CASE ('GFC_FMM');  job_type = GFC_FMM
         CASE ('FE_FMM');   job_type = FE_FMM
         CASE ('QCIM_FMM'); job_type = QCIM_FMM
         CASE DEFAULT; STOP 'fixme'
      END SELECT

      CALL fmm_init_scheme(job_type,LMAX)

   END SUBROUTINE fmm_initial

!-------------------------------------------------------------------------------

   SUBROUTINE fmm_final
      IMPLICIT NONE
      CALL fmm_free_memory
   END SUBROUTINE fmm_final

!-------------------------------------------------------------------------------

   SUBROUTINE fmm_init_qlm_and_V(scheme,dens)

      USE fmm_qlm_builder,     ONLY: fmm_get_raw_qlm
      USE fmm_aux_qlm_builder, ONLY: fmm_get_aux_qlm

      IMPLICIT NONE
      TYPE(fmm_scheme), INTENT(IN) :: scheme
      REAL(REALK),      INTENT(IN) :: dens(:,:)

      NULLIFY (Vff)
      ! Get basic multipole data
      CALL fmm_get_raw_qlm(scheme,dens,LHS_mms,RHS_mms)
      ! Now ensure appropriate prefactors, normalisation and packing
      CALL fmm_get_aux_qlm(scheme,LHS_mms,RHS_mms)
      ! Allocate the far field potential
      CALL fmm_allocate_Vff(scheme)

   END SUBROUTINE fmm_init_qlm_and_V

!-------------------------------------------------------------------------------

   SUBROUTINE fmm_free_memory

      USE fmm_qlm_builder, ONLY: fmm_deallocate_qlm
      IMPLICIT NONE
      IF (ASSOCIATED(Vff)) DEALLOCATE(Vff)
      NULLIFY (Vff)
      CALL fmm_deallocate_qlm(LHS_mms,RHS_mms)

   END SUBROUTINE fmm_free_memory

!-------------------------------------------------------------------------------

   SUBROUTINE fmm_allocate_Vff(scheme)

      IMPLICIT NONE
      TYPE(fmm_scheme), INTENT(IN) :: scheme
      INTEGER(INTK) :: lm_dim, mms_dim, alloc_error

      IF (.NOT.ASSOCIATED(LHS_mms%paras)) STOP 'mms ptrs not set in fmm_main!'
      IF (ASSOCIATED(Vff)) CALL fmm_quit('Vff should NOT be allocated already!')

      mms_dim = SIZE(LHS_mms%paras)
      ! Note we shouldn't use this array for translated, BOXED potentials
      lm_dim = (1+ scheme%raw_LMAX)**2
      IF (scheme%job_type == GFC_FMM) lm_dim = 1

      WRITE(IO1,*) 'Vff: Attempting to allocate',  &
                     MAX(1,lm_dim*mms_dim*8/1000000), 'MB of memory...'
      ALLOCATE (Vff(lm_dim,mms_dim), STAT=alloc_error)
      IF (alloc_error /= 0) WRITE(IO1,*) '... Failed!'

      ! Must zero out since Vff is built additively in general
      Vff(:,:) = zero

   END SUBROUTINE fmm_allocate_Vff

!-------------------------------------------------------------------------------

   SUBROUTINE fmm_get_J_matrix(n_el,dens,J_matrix)

      USE fmm_stats
      USE fmm_scheme_module, ONLY: fmm_get_scheme
      USE fmm_Vff_driver,    ONLY: fmm_get_Vff
      USE fmm_J_builder,     ONLY: fmm_get_J_via_raw_potentials

      IMPLICIT NONE
      CHARACTER(6),  INTENT(IN)  :: n_el
      REAL(REALK),   INTENT(IN)  :: dens(:,:)
      REAL(REALK),   INTENT(OUT) :: J_matrix(:,:)

      TYPE(fmm_scheme), POINTER :: scheme
      CHARACTER(36) :: E
      REAL(REALK)   :: energy, T0, T1

      CALL CPU_TIME(T0)

      CALL fmm_get_scheme(scheme)

      scheme%LHS_mm_range = ELECTRONIC_ONLY
      SELECT CASE (n_el)
         CASE('ONE_EL')
            scheme%RHS_mm_range = NUCLEAR_ONLY
         CASE('TWO_EL')
            scheme%RHS_mm_range = ELECTRONIC_ONLY
         CASE('FULL_J')
            scheme%RHS_mm_range = ALL_MOMENTS
         CASE DEFAULT
            CALL fmm_quit ('require 1, 2, or full J_matrix build!')
      END SELECT
      ! We only have the density on RHS when getting J-matrix
      scheme%LHS_dens = .FALSE.
      scheme%RHS_dens = .TRUE.

      ! Prepare moments and allocate potential
      CALL fmm_init_qlm_and_V(scheme,dens)
      ! Get potential
      CALL fmm_get_Vff(scheme,LHS_mms%paras,RHS_mms,Vff)
      ! Get J-matrix
      CALL fmm_get_J_via_raw_potentials(scheme,Vff,LHS_mms,J_matrix,energy,E)
      WRITE(IO1,'(X,A," = ",E20.12)') E, energy

      CALL fmm_free_memory
      CALL CPU_TIME(T1)
      CALL fmm_TIMTXT('>>> TIME USED in fmm_get_J_matrix', T1-T0, IO1)
      CALL fmm_print_stats

   END SUBROUTINE fmm_get_J_matrix

!-------------------------------------------------------------------------------
! Note that the potential returned by this routine can only be
! contracted with appropriately scaled multipole moments

   SUBROUTINE fmm_get_multipole_potential(mode,npoints,coor,dens,potential)

      USE fmm_stats
      USE fmm_scheme_module, ONLY: fmm_get_scheme
      USE fmm_boundary,      ONLY: fmm_opt_near_field
      USE fmm_Vff_driver,    ONLY: fmm_get_Vff

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: mode
      INTEGER(INTK), INTENT(IN)  :: npoints
      REAL(REALK),   INTENT(IN)  :: coor(3,npoints)
      REAL(REALK),   INTENT(IN)  :: dens(:,:)
      REAL(REALK),   INTENT(OUT) :: potential(:,:)

      TYPE(fmm_scheme), POINTER :: scheme
      REAL(REALK) :: T0, T1
      INTEGER(INTK) :: lmdim

      CALL CPU_TIME(T0)
      CALL fmm_get_scheme(scheme)

      ! Write files for LHS potential grid points
      WRITE(IO1,*) "FMM: npoints = ", npoints
      CALL fmm_write_grid_points(npoints,coor)
      scheme%LHS_mm_range = POTENTIAL_ONLY

      ! We assume files exist for both electronic qlm and nuclear charges
      scheme%RHS_mm_range = ALL_MOMENTS
      !scheme%RHS_mm_range = ELECTRONIC_ONLY
      scheme%LHS_dens = .FALSE.
      scheme%RHS_dens = .TRUE.
      scheme%pack_LHS = .FALSE.

      ! Prepare moments and allocate potential
      CALL fmm_init_qlm_and_V(scheme,dens)

      ! Test if we can skip the near-field interactions
!FIXME
      IF (mode == GFC_FMM) THEN
         CALL fmm_opt_near_field(scheme,LHS_mms%paras,RHS_mms%paras)
      END IF

      ! Get potential
      CALL fmm_get_Vff(scheme,LHS_mms%paras,RHS_mms,Vff)

      ! Note we assume here that the LHS hasn't been rearranged
      lmdim = SIZE(potential,1)
      IF (SIZE(potential,2) /= SIZE(Vff,2)) CALL fmm_quit('bounds: potential')
      potential(:,:) = Vff(1:lmdim,:)

      CALL fmm_free_memory
      CALL CPU_TIME(T1)
      CALL fmm_TIMTXT('>>> TIME USED in fmm_get_potential', T1-T0, IO1)
      CALL fmm_print_stats

   END SUBROUTINE fmm_get_multipole_potential

!-------------------------------------------------------------------------------

   SUBROUTINE fmm_write_grid_points(npoints,coor)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN) :: npoints
      REAL(REALK),   INTENT(IN) :: coor(3,npoints)

      INTEGER(INTK) :: i

      ! Write grid points to disk
      OPEN(UNIT=IO2, FILE='grid.fmm', STATUS='REPLACE',  &
           ACCESS='SEQUENTIAL', FORM='UNFORMATTED')
      REWIND(IO2)
      WRITE(IO2) npoints
      DO i = 1, npoints
         WRITE(IO2) coor(1,i), coor(2,i), coor(3,i)
      END DO
      CLOSE(UNIT=IO2, STATUS='KEEP')

   END SUBROUTINE fmm_write_grid_points

!-------------------------------------------------------------------------------

END MODULE fmm_main

